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FULLY  MULTIDIMENSIONAL  FLUX-CORRECTED  TRANSPORT 


I.  INTRODUCTION:  FCT  DEFINED 

Consider  the  following  system  of  equations 

w,  + fx  — 0 ( 1 ) 

where  w and  / are  vector  functions  of  independent  variables  xand  r.  A simple  example  of  such 

a system  of  equations  would  be  the  one  dimensional  equations  of  ideal,  inviscid  fluid  flow: 


w = 

p 

pv 

; / = 

pv 

pv2  -1-  P 

pF 

pv£  + P v 

where  p,  v,  P and  £are  the  fluid  density,  velocity,  hydrostatic  pressure  and  specific  total  energy 
respectively. 

We  shall  say  that  a finite  difference  approximation  to  Eq.  (1)  is  in  conservation  (or  "flux") 
form  when  it  can  be  written  in  the  form 

<+1  = A*,-1  [^+(1/2)  - f,— <i/2)]  (2) 

Here  w and  / are  defined  at  the  spatial  grid  points  x,  and  temporal  grid  points  t",  and  A.v,  = 

y (x,+1  — x,^).  The  /r,+(l/2)  are  called  transportive  fluxes,  and  are  functions  of  / at  one  or 

more  of  the  time  levels  t".  The  functional  dependence  of  F on  f defines  the  integration  scheme 
(leapfrog,  Lax-Wendroff,  Crank  Nicholson,  donor  cell,  etc.). 

It  is  well  known  that  higher  order  (order  2 and  above)  schemes  for  numerically  integrat- 
ing Eq.  (1)  suffer  from  dispersive  "ripples"  in  w,  particularly  near  steep  gradients  in  w.  Lower 
order  schemes,  such  as  donor  cell,  Lax-Friedrichs,  or  high  order  schemes  with  a zeroth  order 
diffusion  added,  produce  no  ripples  but  suffer  from  excessive  numerical  diffusion.  Flux- 
corrected  transport  (FCT)  is  a technique  developed  by  Boris  and  Book  [1-3]  which  embodies 
the  best  of  both  of  the  above  worlds.  In  its  simplest  terms,  FCT  constructs  the  net  transportive 
flux  point  by  point  ( non  linearly)  as  a weighted  average  of  a flux  computed  by  a low  order 
scheme  and  a flux  computed  by  a high  order  scheme.  The  weighting  is  done  in  a manner  which 
insures  that  the  high  order  flux  is  used  to  the  greatest  extent  possible  without  introducing  rip- 
’Manuscript  submitted  May  17,  1978 
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pies  (overshoots  and  undershoots).  This  weighing  procedure  is  referred  to  as  "flux-correction" 
or  "flux-limiting"  for  reasons  which  shall  become  clear  later.  The  result  is  a family  of  transport 
algorithms  capable  of  resolving  moving  contact  discontinuities  over  3-4  grid  points,  and  shock 
fronts  over  2 grid  points,  without  undershoot  or  overshoot  [1-3).  Formally,  the  procedure  is  as 
follows: 

1)  Compute  F,+(1/2),  the  transportive  flux  given  by  some  low  order  scheme  guaranteed  to 
give  monotonic  (ripple-free)  results  for  the  problem  at  hand 

2)  Compute  F,+(1/2),  the  transportive  flux  given  by  some  high  order  scheme 

3)  Define  the  "antidiffusive  flux": 

A,+t\/2)  = ^,  + (1/2)  ~ ^ + (1/2) 

4)  Compute  the  updated  low  order  ("transported  and  diffused")  solution: 

*,;</=  wn_  aX(-i  | Ff+um  - /r,-(i/2)] 

5)  Limit  the  /4,+(1/2)  in  a manner  such  that  w',  + l as  computed  in  step  6 below  is  free  of 
overshoots  and  undershoots: 

^<+(1/2)  = C<+(l/2)  ^,+(I/2)<  0 ^ C+(l/2)  ^ 1 

6)  Apply  the  limited  antidiffusive  fluxes: 

<+l  = w!d~  A*,’1  [>f+<i/2>  ~ ^.-(1/2)] 

The  critical  step  in  the  above  is,  of  course,  step  5 which  will  be  discussed  shortly.  In  the 
absence  of  the  flux  limiting  step  (-4, +0/2)  = <4,+<i/2)).  w,"  + l would  simply  be  the  time-advanced, 
high  order  solution. 

We  note  that  this  definition  of  FCT  is  considerably  more  general  than  has  been  given  pre- 
viously. 

II.  MULTIDIMENSIONAL  FLUX-CORRECTED  TRANSPORT 

Before  proceeding  to  a discussion  of  flux  limiting,  let  us  see  how  the  procedure  given 
above  might  be  implemented  in  multidimensions.  An  obvious  choice  would  be  to  use  a 


Strang-type  time-splitting  procedure  [4]  when  it  can  be  shown  that  the  equations  allow  such  a 
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technique  to  be  used  without  serious  error.  Indeed,  such  a procedure  may  even  be  preferable 
from  programming  and  time-step  considerations.  However,  there  are  many  problems  for  which 
time-splitting  produces  unacceptable  numerical  results,  among  which  are  those  involving 
incompressible  or  nearly  incompressible  flow  fields.  This  technique  is  straightforward  and  shall 
not  be  discussed  here.  Instead  we  consider  as  an  example  the  fully  two-dimensional  set  of  equa- 
tions 


H>,  + fx  + gy  “ 0 (3) 

where  w,  /,  and  gare  vector  functions  of  x,  y , and  t.  In  finite  difference  flux  form  we  have 

w"'j ' = w"j  - A y~}  |^,+(i/2).y  - Fi-am.j  + gi.j+(  1/2)  ~ Gi.  v-d/2)]  (4) 

where  now  w,  / and  g are  defined  on  spatial  grid  points  x„  y,  at  time  levels  t",  and  A V,  j is  a 

two  dimensional  area  element  centered  on  grid  point  (/,  j).  Now  there  are  two  sets  of  tran- 
sport! ve  fluxes  F and  G , and  the  FCT  algorithm  proceeds  as  before: 

1)  Compute  F,l+(i/2),  j and  G,f->+(1/2)  by  a low  order  monotonic  scheme 

2)  Compute  F,+a/2). , and  G,Hj+( i/2)  by  a high  order  scheme 

3)  Define  the  antidiflfusive  fluxes: 

^<+(1/2).;  — ^/+(i/2). ; “ Fh-nm.j 

Ai,j+ d/2)  = Gu+un)  ~ Gtj+um 

4)  Compute  the  low  order  time  advanced  solution: 

w!fj  “ w"j  - A V Cj  [^(1/2).,  “ ^-(1/2),  J + G,f7  + (1/2)  - G,fy_(1/2)] 

5)  Limit  the  antidiflfusive  fluxes 

Af+(\ll),j  = -4,  + (1/2),  j C+O/2  ).J  0 ^ C + (l/2  ).J  ^ 1 

dfj+U/2)  = ^i.J+W2).J  Q j+(\/2)  0 < Ctij+ (1/2)  < 1 

6)  Apply  the  limited  antidiflfusive  fluxes: 

w"j'  = wu  - AV“J  " Ah-(\/2).j  + ^Fj+d/ 2)  ~ /5y— d/2)] 

As  can  be  easily  seen,  implementation  of  FCT  in  multidimensions  is  straightforward  with 
the  exception  of  Step  5,  an  algorithm  for  which  is  the  subject  of  this  paper.  First,  however,  let 
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us  see  how  flux  limiting  is  presently  implemented  in  one  spatial  dimension. 


III.  FLUX  LIMITING  IN  ONE  SPATIAL  DIMENSION  - 
THE  ORIGINAL  ALGORITHM 

The  original  algorithm  for  flux-limiting  in  one  dimension  was  given  by  Boris  and  Book 
[1].  In  our  notation  it  is: 


<+(1/2)  — ^(+<1/2)  max 


0,  min  k+(1/2)| 


5'+<l/2)  ( w'+2  - w!+]\  A.V,  + 1,  A,+(l/2)  ( W,"'  - w't 1 1 A.V,_iJ 


(5) 


where 


o _ + 1 if  /U+(i/2)  ^ 0 
(+0.2)  _ j jf  Al+(]/2)  < o 

The  intent  of  this  formula  is  that  antidiffusive  fluxes  should  neither  create  new  extrema,  nor 
accentuate  already  existing  extrema,  in  the  transported  and  diffused  solution  w'd.  That  the 
above  formula  does,  in  fact,  perform  precisely  this  task  can  be  verified  by  the  reader  with  rela- 
tive ease.  We  shall  examine  here  some  of  the  less  obvious  properties  of  this  very  powerful,  yet 
simple,  formula.  In  the  process  we  shall  gain  insight  into  which  of  these  properties  we  shall 
wish  to  carry  over  into  a multidimensional  flux  limiter.  We  first  observe  that  certain  quantities 
do  not  appear  in  the  above  formula:  1)  w,"^  - w,"',  the  first  difference  of  w"'at  the  point  where 
the  antidiffusive  flux  -4,+<l/2)  is  evaluated;  and  2)  antidiffusive  fluxes  other  than  A,+( l/2).  This 
last  property  is  the  most  notable  since  there  are  conceivably  two  fluxes  directed  into  or  out  of  a 
cell.  A formula  guaranteeing  that  the  two  fluxes  acting  in  concert  shall  not  create  ripples  would 
apparently  require  a knowledge  of  both.  We  shall  return  to  this  point  momentarily. 

In  Figure  1 we  show  the  eight  possible  configurations  of  w"1  in  the  neighborhood  of  a 
positive  /f,+d/2)  ( directed  to  the  right  in  our  diagrams).  Configurations  1-4  show  the  "normal'1 
situation,  with  A,+^n)  having  the  same  sign  as  w,"',  - w,"'  (as  might  be  expected  of  an 
antidiffusive  flux  ).  We  note  that  if  cither  w''l2  — w,  + | or  w,"1  — wjt\  has  a sign  opposed  to  that 
^(+0/2)'  as  in  configurations  2-4,  the  antidiffusive  flux  ^4 , + r i /2>  is  completely  canceled.  This, 


4 


however,  is  in  total  agreement  with  the  stated  intent  of  Eq.  (5)  since  otherwise  configuration  2 
would  allow  accentuation  of  an  existing  maximum,  configuration  3 accentuation  of  an  existing 
minimum,  and  configuration  4 accentuation  of  both.  In  the  remaining  configuration  1,  the  flux 
limiter  (5)  will  reduce  the  magnitude  of  /4(+(1/ 2)  sufficie''"-'  to  guarantee  that  neither  a max- 
imum at  grid  point  i + 1 nor  a minimum  at  grid  point  / will  be  formed,  again  in  precise  agree- 
ment with  its  stated  intent. 

Configurations  5-8  are  identical  to  configurations  1-4,  respectively,  except  that  sign  of 
— wld  has  been  reversed  (The  "antidiffusive  fluxes"  are  now  directed  down  the  gradient  in  w'J). 
Since  the  sign  of  w,+,  - w'd  does  not  enter  into  the  flux  correction  formula  (5),  the  results  of 
the  formula  are  identical  to  those  for  the  previous  four  cases:  the  antidiffusive  fluxes  are  can- 
celed for  configurations  6-8  and  limited  in  configuration  5 to  the  extent  necessary  to  prevent  a 
new  maximum  at  grid  point  / + 1 and  a new  minimum  at  grid  point  /.  Examination  of 
configurations  6-8  reveals  that  -4,  + <i/2>  actually  presented  no  hazard  insofar  as  extrema  creation 
or  enhancement  (at  least  in  moderation).  Certainly  there  was  no  cause  for  completely  cancel- 
ing the  flux.  Even  in  configuration  5 the  flux  may  have  been  limited  to  a greater  extent  than 
necessary.  At  first  it  would  seem  that  configurations  5-8  represent  errors  introduced  by  the 
simplicity  of  the  flux  limiting  formula  (5).  However,  extensive  tests  by  this  author  indicate 
that  in  the  relatively  rare  instances  in  which  configurations  5-8  occur  in  practice,  the  "errors" 
introduced  by  Eq.  (5)  represent,  in  fact,  the  correct  action  to  take  in  terms  of  producing  accu- 
rate profiles  in  w"+l.  More  importantly,  they  represent  the  mechanism  by  which  Lq.  (5)  can 
guarantee  that  ripples  are  not  formed  under  any  circumstances,  as  we  shall  see  presently. 

Consider  two  antidiffusive  fluxes,  acting  in  concert,  attempting  to  produce  or  accentuate 
an  extremum.  We  therefore  have  -4,+()/2)  and  -4,_(|/2)  either  both  directed  toward,  or  both 
directed  away  from  grid  point  i.  We  see  from  Figure  1 that,  in  general,  an  antidiffusive  flux 
directed  opposite  to  the  gradient  in  H>'rfwill  be  completely  canceled.  Therefore  the  only  cases  of 
fluxes  acting  in  concert  that  we  need  be  concerned  with  are  those  where  two  adjacent  fluxes  are 
both  parallel  to  the  local  gradients  in  w'd.  These  are  precisely  the  cases  of  already  existing 


extrema,  in  which  case  both  fluxes  will  be  canceled  (as  in  configurations  2-4).  This  is  the  reason 
that  Eq.  (5)  needs  no  information  on  any  antidiffusive  flux  other  than  ^ +< i/2>. 

In  Figure  2 we  see  that  the  above-mentioned  assumptions  regarding  antidiffusive  fluxes 
acting  in  concert  break  down  completely  in  multidimensions.  It  is  possible  in  more  than  one 
dimension  for  more  than  one  antidiffusive  flux  to  be  directed  into  or  out  of  a cell,  all  of  these 
fluxes  being  directed  parallel  to  the  local  gradient  in  w'd,  without  that  cell  being  an  already 
existing  extremum.  Therefore  the  problem  of  dealing  with  multiple  antidiffusive  fluxes  acting  in 
concert  cannot  be  avoided  by  simply  canceling  all  fluxes  antiparallel  to  the  local  gradient  in  w'd. 
It  is  clear  then  that  any  formula  which  purports  to  perform  flux  limiting  in  more  than  one 
dimension  without  resort  to  time  splitting  must  contain  information  about  antidiffusive  fluxes 
other  than  the  one  being  limited. 

IV.  FLUX  LIMITING  IN  ONE  SPATIAL  DIMENSION  - 
AN  ALTERNATIVE  ALGORITHM 

We  describe  here  in  one  spatial  dimension  an  alternative  flux  limiting  algorithm  which 
generalizes  easily  to  multidimensions  and  which,  even  in  one  dimension,  exhibits  a superiority 
to  the  limiter  described  in  the  previous  section  (Eq,  (5))  with  regard  to  peaked  profiles. 

Referring  to  Figure  3,  we  seek  to  limit  the  antidiffusive  flux  -4,+(i/2>  such  that 

^i  + d/2)  = C + (I/2)  /L  + U/2>>  0 ^ C + tl/2)  ^ 1 (6) 

and  such  that  /L+d/2)  acting  in  concert  with  will  not  allow 

w"*'  = wtd  — A.v,-1  |/4,c+(|/2)  - <4,-d/:)j 

to  exceed  some  maximum  value  w,max  nor  fall  below  some  minimum  value  w,m,n.  We  leave  the 
determination  of  w,ma’1  and  w,mm  until  later. 

We  define  three  quantities: 

P,+  = the  sum  of  all  antidiffusive  fluxes  into  grid  point  / 


= max  (0,  - min  (0,  -4,+<i/2>) 

(7) 

Q*  = ( H',max  — H’/'O  \X, 

(8) 

mind,  Q,+/P,+)  if  P,+  > 0 

Ki  o if  p;  = o 

(9) 
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Assuming  that  wm"  > w,'d  (it  must  be),  all  three  of  the  above  quantities  are  positive  and  R f 
represents  the  least  upper  bound  on  the  fraction  which  must  multiply  all  antidiffusive  fluxes  into 
grid  point  i to  guarantee  no  overshoot  at  grid  point  /. 


Similarly  we  define  three  corresponding  quantities: 

p - = the  sum  of  all  antidiffusive  fluxes  away  from  grid  point  i 

= max  (0,  A(+,i/2))  - min  (0,  <4,-U/2)) 

Q- - (w’d  - w,m'n)  Ax, 


R, 


min(l,  Qr/PD  if  P,~  > 0 
0 if  Pf  = 0 


(10) 

(11) 

(12) 


Again  assuming  that  w,min  < w*  we  find  that  R,~  represents  the  least  upper  bound  on  the 
fraction  which  must  multiply  all  antidiffusive  fluxes  away  from  grid  point  i to  guarantee  no  un- 
dershoot at  grid  point  /. 

Finally  we  observe  that  all  antidiffusive  fluxes  are  directed  away  from  one  grid  point  and 
into  an  adjacent  one.  Limiting  will  therefore  take  place  with  respect  to  undershoots  for  the 
former  and  with  respect  to  overshoots  for  the  latter.  A guarantee  that  neither  event  comes  to 


pass  demands  our  taking  a minimum: 

(min  {Rj+i,  R~)  *f  'L+a m ^ ® ^3) 

Cj+u/2)  _ |mjn  (/?  +,  R~+x)  if  A, +0/2)  < 0 

Furthermore,  we  shall  call  upon  our  previously  described  experience  with  the  original  flux 
limiter  and  set 

A, +<i/2)  = 0 if  A, +(i/2)  (wHy  - w<*)  < 0 (H) 

and  either  4,+(i/2)(wi+2  ~ w!i0  < ® 

or  -4, +(i/2)  (w/rf—  <0  <0 

In  practice  the  effect  of  Eq.  (14)  is  minimal  and  is  primarily  cosmetic  in  nature.  This  is 
because  cases  of  antidiffusive  fluxes  directed  down  gradients  in  w,d  are  rare,  and  even  when 
they  occur  usually  involve  flux  magnitudes  that  are  small  compared  to  adjacent  fluxes.  If  Eq.  (14) 
is  used,  it  should  be  applied  before  Eq.  (6)  through  (13). 
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We  come  now  to  a determination  of  the  quantities  tv,max  and  tv  mi"  in  Eqs.  (8)  and  (11). 


A safe  choice  is 

tv,max  = max  (vt’/'i,,  wjd,  tv,+ , ) (15) 

tv,min  = min  ( w'd  | , w'J,  vv/^i ) (16) 

This  choice  will  produce  results  identical  with  those  of  Eq.  (5)  in  one  dimension,  including  the 

occurrence  of  the  "clipping"  phenomenon  to  be  mentioned  shortly. 

A better  choice  is: 


tv®  = max  ( tv",  tv,"4) 


tv,max  = max  ( w,a_ , , 

<.  <+i ) 

(17) 

w,h  = min  ( tv,” 

tv,"4) 

H’,mm  = min  (w,4.]. 

tv,4,  tv4.,) 

(18) 

This  choice  allows  us  to  look  back  to  the  previous  time  step  for  upper  and  lower  bounds  on 
tv/,+1. 

It  is  clear  that  these  two  methods  of  determining  tv,max  and  tv,mm  represent  only  a small 
subset  of  possible  methods.  The  alternative  flux  limiter  described  in  equations  (6)  through  (14) 
admits  of  any  physically  motivated  upper  and  lower  bound  on  tv,"  1 supplied  by  the  user,  intro- 
ducing a flexibility  unavailable  with  the  original  flux  limiter  (5).  However,  with  the  exception 
of  one  example  in  the  next  section  (which  shows  graphically  the  potential  power  of  this  flexibil- 
ity), we  shall  henceforth  use  Eq.  (17)  and  (18)  to  evaluate  tvmin  and  tvmax  in  one  dimension. 

V.  COMPUTATIONAL  EXAMPLES  - ONE  DIMENSION 

We  consider  one  dimensional  passive  convection  in  a constant  velocity  field.  We  have  Eq. 
(1)  with  tv  = p and  / = pv  with  v = constant.  We  choose  our  transport  algorithm  to  be  that 
given  in  [3]  for  LPE  Shasta.  On  the  standard  square  wave  tests  we  find  that  our  results  for  the 
original  flux  limiter  (5)  and  for  the  alternative  flux  limiter  (6)  through  (14)  are  identical  to 
within  round-off  (the  same  is  true  for  traveling  shock  waves  in  the  coupled  one  dimensional 
equations  of  ideal  inviscid  fluid  flow).  To  find  differences  between  the  limiters  in  one  dimen- 


sion we  must  look  to  passive  convection  of  peaked  profiles.  We  choose  the  problem  given  by 

Forester  [5],  a guassian  of  half-width  2A.v.  In  Figure  4 we  show  the  results  after  600  iterations 
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for  the  trivial  case  v = 0.  On  the  left  we  see  the  fimilar  "clipping"  phenomenon  with  the  origi- 
nal flux  limiter,  caused  by  a zeroth  order  diffusion  term  in  the  low  order  portion  of  the  LPE 
Shasta  algorithm.  This  diffusion  term  causes  the  peak  in  w'^to  be  smaller  than  the  peak  in  w", 
leaving  the  original  flux  limiter  (5)  with  no  way  of  resurrecting  the  original  peak.  This  process 
occurs  repeatedly,  eventually  leaving  the  characteristic  three  point  top.  The  alternative  flux 
limiter,  shown  on  the  right,  "remembers"  the  old  value  of  the  peak  and  is  able  to  resurrect  it 
each  time  step. 

In  Figure  5 we  show  the  same  problem  after  600  iterations  but  this  time  for  e = v Ar/Ax 
= 0.1.  We  see  that  clipping  occurs  with  both  flux  limiters,  but  to  a lesser  extent  with  the  alter- 
native flux  limiter  (6)  through  (14). 

At  this  point  we  removed  the  flux  limiter  entirely  and  again  ran  the  problem  600  itera- 
tions with  e = 0.1.  The  results  convinced  us  that  the  amplitude  and  phase  properties  of  the 
high  order  portion  of  LPE  Shasta  were  incapable  of  resolving  the  high  wave  numbers  of  which 
the  gaussian  is  composed.  Consequently  it  was  decided  to  switch  to  a higher  order  algorithm,  a 
leapfrog-trapezoidal  transport  algorithm  which  uses  eighth  order  spatial  differences.  The  algo- 
rithm is,  then,  second  order  accurate  in  time  and  eighth  order  in  space,  with  an  amplification 
factor  that  is  effectively  unity  across  the  entire  Fourier  spectrum,  and  phase  properties  consider- 
ably better  than  those  of  fourth  order  algorithms.  The  leapfrog  portion  of  this  algorithm  is 
identical  to  the  eighth  order  Kreiss-Oliger  scheme  [6].  A fourth  order  version  of  this  same 
algorithm  was  used  later  in  the  two-dimensional  solid  body  rotation  tests.  We  ran  the  gaussian 
test  problem  again  600  iterations  with  e = 0.1  with  no  flux  limiter  and  were  convinced  that  the 
algorithm  did  indeed  have  the  resolving  power  necessary  to  do  the  problem.  A low  order 

scheme,  donor  cell  plus  a zeroth  order  diffusion  term  with  coefficient  y,  was  added  to  com- 

O 

plete  the  FCT  algorithm,  which  we  dub  2-8  leapfrog-trapezoidal.  A more  detailed  description 
of  this  algorithm  is  found  in  the  appendix. 

Figure  6 shows  the  results  of  2-8  leapfrog  trapezoidal  run  600  iterations  with  e = 0.1  with 

both  the  original  and  alternative  flux  limiters.  The  results  are  better  than  those  in  Figure  5, 
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and  again  the  alternative  flux  limiter  proves  superior,  but  nontheless  disappointing.  The  clip- 
ping would  appear  to  be  due  entirely  to  the  flux  limiters,  not  to  the  phase  or  amplitude  proper- 
ties of  the  high  order  scheme. 

A careful  examination  of  exactly  what  happens  to  a one  point  peak  in  a finite  difference 
code  reveals  the  real  source  of  the  above  problem.  Consider  a profile  with  a local  peak  at  grid 
point  i in  passive  convection  at  constant  velocity  > 0.  At  each  succeeding  lime  step  the  func- 
tion value  at  grid  point  / will  decrease  and  that  at  / + 1 will  increase  (Figure  7).  Eventually 
they  will  both  reach  some  intermediate  value,  and  the  actual  original  peak  value  will  not  appear 
anywhere  on  the  grid,  since  it’s  position  now  lies  midway  between  two  grid  points.  At  this 
point  even  the  new  flux  limiter  (6)  through  (14),  (17)  and  (18),  has  lost  the  information  it 
needs  to  allow  the  peak  to  be  resurrected  in  suceeding  time  steps,  and  will  "clip"  the  new  peak 
at  grid  point  / + 1 as  it  tries  to  form,  based  on  the  assumption  that  it  is,  in  fact,  an  overshoot. 
The  effect  is  magnified,  since  the  clipping  itself  introduces  phase  errors  in  succeeding  iterations, 
the  net  result  being  the  profiles  dipicted  in  Figure  6. 

It  is  clear,  then,  that  if  we  are  to  successfully  treat  a one-point  extremum  within  the  con- 
text of  FCT  we  must  use  information  other  than  just  the  grid  point  values  themselves.  In  what 
follows  we  shall  utilize  the  flexibility  of  the  alternative  flux  limiter  to  use  as  w,max  and  w,mm  any 
values  that  we  choose.  In  Figure  8 we  show  a possible  way  of  extracting  information  about 
extrema  which  do  not  lie  exactly  on  a grid  point  at  the  time.  Basically  we  define  to  be 

the  w value  at  the  intersection  of  the  line  segments  formed  by  connecting  the  point  (x,_|,  w/l) ) 
with  (x„  w'd)  and  the  point  (x,+1,  *>/£,)  with  (x,+2.  w,+2).  If  the  x coordinate  of  this  intersec- 
tion lies  between  x,  and  x,  + 1,  then  we  consider  this  to  be  an  allowable  wmax  or  wm,n  for 

either  w"  + 1 or  w,"+V-  We  now  have 

w“  = max  ( w",  w,"/) 

w,max  = max  w“,  w,°+l,  tv,^f,k/2),  w,pe(aj‘/2) ) (19) 

w,h  = min  (w",  tv,"0 

wmin  = min  «„  w!>,  w/V„  w,?tf/2) ) (20) 

Equations  (19)  and  (20)  together  with  equations  (6)  through  (14)  now  determine  the 
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alternative  flux  limiter  (for  this  section  only). 

In  Figure  9 we  show  the  results  of  using  Equations  (19)  and  (20)  to  determine  w,max  and 
w,mm  on  the  gaussian  test  problem  run  600  iterations  with  e = 0.1.  Clearly  the  problem  has 
been  solved  — we  recover  the  gaussian  profile  with  no  dispersive  ripples  and  minimal  loss  in 
amplitude.  We  have  not  performed  this  test  merely  to  show  the  power  of  the  extrapolation 
technique  just  described  to  determine  w,max  and  w,min.  Rather  this  calculation  serves  to  show  the 
power  of  using  information  other  than  that  available  on  the  one  dimensionht  grid.  In  multidimen- 
sional flux  limiting,  this  information  comes  from  the  other  coordinate  directions,  as  we  shall 
see. 

VI.  FLUX  LIMITING  IN  MULTIDIMENSIONS 

The  alternative  flux  limiting  algorithm  presented  in  section  IV  generalizes  trivially  to  any 
number  of  dimensions.  For  the  sake  of  completeness  we  present  here  the  algorithm  for  two 
spatial  dimensions. 

Referring  to  Figure  10,  we  seek  to  limit  the  antidiffusive  fluxes  /4,+(1/2 >,  / and  A, , ;+(|/2) 
such  that 

-4  . + (1/2).;  = C,+(l/2),y  ‘4-  + (l/2),7  0 < C,  + (l/2) .j  ^ 1 
AU+W2)  = C.j  + O/2)  Ai.j+ (1/2)  0 < C;,y+(  1/2)  ^ 1 

and  such  that  Af+nn)^,  AfL(U2).j,  Afj+( 1/2),  and  A^_iU2)  acting  in  concert  shall  not  cause 

O1  = w!dj  - & K7j  [^,  + (1/2).,  “ ^--(1/2).;  + ^ w + (l/2)  “ ^/C>-(l/2)] 

to  exceed  some  maximum  value  w,max  nor  fall  below  some  minimum  value  w,mjn. 

Again  we  compute  six  quantities  completely  analogous  to  those  computed  in  Eq.  (7) 
through  (12): 

P+j  = the  sum  of  all  antidiffusive  fluxes  into  grid  point  (/j  j)  (7  ) 

- max  (0,  ^,-d/2).  y)  ~ min  '0.  A,+(V2).j) 

+ max  (0,  /l(iy_,i/2))  - min  '0,  -4,,y+(i/2)) 

q+j  = ( - w?j)  A V,  j (8') 

|min(I,  Q+j/P+j  if  P*j  > 0 ( ) 

X'--|  0 lf^-0 
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Pjj  “ the  sum  of  all  antidilTusive  fluxes  away  from  grid  point  (/,  J) 


- max  (0,  Ai+im),j)  - min  (0,  v4(_<i/2).y) 
+ max  (0,  Au+im))  - min  (0,  AU_W1)) 
QCj-("u~  A y,  j 

|min(l,  QCj/P-j)  if  P-j  > 0) 

' | 0 if^-y-0) 

Equation  (13)  becomes 


C/+(i/2),y 
Ci.  y+a/2) 

while  Eq.  (14)  becomes 


(min  (RjXi  j,  /?,y)  if  ■^/+(i/2).y  ^ 0 
|min  (R+j,  Rj+i,j)  if  ^+0/2).  y < 0 

(min  (R+j+ 1,  /?,";)  if  -4,.  7+0/2)  ^ 0 
(min  (7?,+y,  R,~j+i ) if  y-Ki/2)  < 0 


(10') 

on 

(12') 


(13') 


^/+(i/2),y  “ 0 if  -4,+(i/2).y  (w/+i.y  “ <y)  < 0 
and  either  Al+(l/1)ij  (w/+2.y  ~ w/£i.y)  < 0 
ot  ^1+0/2),  y «y  ~ w'-i.y)  < 0 
^-.y+(i/2)  “ 0 if  Au+(U2)  (h - w,*7)  < 0 
and  either  /f,,y+(1/2)  (<,+2  - <% ) < 0 

or  au+u/2)  (»u  ~ <y-i>  < 0 04') 

and  Eq.  (17)  and  (18)  become 


w°j  - max  (*,",,  w'fj) 

<7*  - max  (<-l.y,  W°j.  <y_,.  < y+l  ] 

w,fy  - min(w,^,  wjfj) 

<yin  - min(u’,A_1J,  w£.1>y,  w,?y_,,  w,?y+1) 

Again,  the  effect  of  Eq.  (14')  is  minimal,  but  if  it  is  used  it  should  be  applied  before  Eq. 
(6')  through  (13').  Note  that  our  search  for  w,™fx  and  wfin  now  extends  over  both  coordinate 
directions.  Where  finite  gradients  exist  in  both  directions,  this  procedure  will  allow  us  to  stop 
the  clipping  phenomenon  in  regions  where  a peak  exists  with  respect  to  one  coordinate  direc- 
tion but  not  in  the  other,  as  we  shall  see  in  the  next  section. 


(17') 

(18') 


VII.  COMPUTATIONAL  EXAMPLES  - TWO  DIMENSIONS 

We  choose  as  our  two  dimensional  test  problem  that  of  solid  body  rotation.  That  is,  we 
have  Eq.  (3)  with  / - w vx,  g - w v,,  v,  - -fl(y  - y0),  and  v,  - Cl(x  - x„).  Here  fl  is  the 
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(constant)  angular  velocity  in  radians/sec  and  (x„,  y„)  is  the  axis  of  rotation.  The  configuration 
is  shown  in  Figure  11.  The  computational  grid  is  100  x 100  cells.  Ax  = Ay,  with  counterclock- 
wise rotation  taking  place  about  grid  point  (50,  50).  Centered  at  grid  point  (50,  75)  is  a cylinder 
of  radius  15  grid  points,  through  which  a slot  has  been  cut  of  width  5 grid  points.  The  lime 
step  and  rotational  speed  are  chosen  such  that  628  time  steps  will  effect  one  complete  revolu- 
tion of  the  cylinder  about  the  central  point.  A perspective  view  of  the  initial  conditions  is 
shown  in  Figure  12. 

Our  high  order  scheme  for  the  following  tests  is  a fully  two  dimensional,  fourth  order  in 
space,  second  order  in  time  leapfrog-trapezoidal  scheme,  the  leapfrog  step  of  which  is  a two 
dimensional  fourth  order  Kreiss-Oliger  scheme  [61.  The  low  order  scheme  is  simply  two 
dimensional  donor  cell  plus  a two  dimensional  zeroth  order  diffusion  term  with  diffusion 


coefficient 


A more  detailed  description  of  this  algorithm  is  found  in  the  appendix. 


We  wish  to  emphasize  that  the  only  difference  between  calculations  in  the  following  com- 
parisons is  in  the  flux  limiting  stage  itself.  The  high  order  fluxes,  low  order  fluxes,  and  hence 
the  (unlimited)  antidiffusive  fluxes  are  all  computed  in  the  full  two  dimensions  without  using 
time-splitting.  In  each  case  we  are  comparing  the  fully  two  dimensional  flux  limiter  given  by 
Eq  (6')  through  (14'),  (17')  and  (18  ) with  a time  split  application  of  Eq.  (5).  Time  splitting  is 
the  only  way  that  Eq.  (5)  may  be  utilized  in  a multidimensional  problem.  Note  that  in  the 
latter  case  we  are  not  time  splitting  the  entire  transport  operator,  but  only  the  flux  limiter  (5). 
In  this  way  we  are  testing  only  the  limiters  themselves. 

In  Figure  13  we  show  a perspective  view  of  the  two  calculations  after  revolution  (157 


iterations).  Figure  14  presents  a comparison  of  the  results  of  the  two  calculations  for  one  full 
revolution  (628  cycles).  Two  features  are  obvious.  The  first  is  a much  greater  filling-in  of  the 
jiot  with  the  lime  split  Eq.  (5)  than  with  the  fully  two  dimensional  flux  limiter.  The  second  is 
the  loss  of  the  bridge  connecting  the  two  halves  of  the  cylinder  in  the  case  of  the  time-split 
application  of  Eq.  (5).  Less  obvious  is  the  lack  of  dipping  of  the  peaked  profiles  defining  the 
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front  surface  of  the  cylinder  for  the  case  of  the  fully  multidimensional  limiter.  Clearly  this  is 
due  to  the  fact  that  the  multidimensional  flux  limiter  can  look  in  both  directions  to  determine 
whether  or  not  a genuine  maximum  exists.  Note  that  there  are  two  factors  working  in  favor  of 
the  fully  multidimensional  flux  limiter:  1)  the  ability  to  look  in  both  directions  to  find  minima 
and  maxima,  as  just  mentioned;  and  2)  the  ability  to  scan  both  w,"y  and  wlJj  to  find  maxima  and 
minima.  Both  of  these  factors  are  responsible  for  the  improved  profiles. 

VIII.  THE  STRIATIONS  CODE  - A TWO  DIMENSIONAL  INCOMPRESSIBLE 
FLUID  CODE  USING  FULLY  MULTIDIMENSIONAL  FCT 

A two  dimensional  ( x - y)  plasma  cloud  initialized  in  a region  of  constant  magnetic  field 
B„  directed  along  the  z axis,  with  an  externally  imposed  electric  field  E0  directed  along  the  x 
axis  will  tend  to  drift  in  the  E0  x B0  direction  (along  the  negative  y axis)  (see  Figure  15).  If 
the  ion-neutral  collision  frequency  is  finite,  Pedersen  conductivity  effects  will  produce  polariza- 
tion fields  which  tend  to  shield  the  inner  (more  dense)  regions  of  the  cloud  from  E„,  causing 
this  inner  portion  of  the  cloud  to  drift  more  slowly  than  the  outer  portions  of  the  cloud.  This 
results  in  a steepening  of  gradients  on  the  back  side  of  the  cloud.  Arguments  similar  to  those 
above,  applied  to  infinitesimal  perturbations  imposed  upon  this  back  side  gradient,  show  that 
the  back  side  of  the  cloud  is  physically  unstable  to  perturbations  along  x.  For  a detailed 
description  of  this  problem,  see  [7], 

The  equations  of  motion  for  the  electron  fluid  are: 

(bNJbt)  + VA(JVf  V,)  = 0 (21) 

VA  • (NeV^f)  = E„  • VjN,  (22) 

V,  — — (c/B0)  x z*  (23) 

Here  Ne,  V,,  and  'P  are  the  electron  density,  electron  velocity  and  perturbation  electric 

d d 

field  potential  respectively,  and  V±  is  the  two  dimensional  divergence  operator  x — + y — . 

The  magnitudes  of  B„  and  E0  are  0.5  gauss  and  5 millivolts  per  meter  respectively.  Our  rest 
frame  here  is  that  of  the  (c/B0) E„  x z velocity.  A few  trivial  vector  identities  will  convince  the 
reader  that  VAVr  - 0,  meaning  that  the  electrons  move  incompressibly.  Clearly  time-splitting 
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the  transport  operator  would  be  disastrous  here,  and  a fully  two  dimensional  scheme  is  re- 
quired. 

Eq.  (22)  is  solved  for  ^ using  an  elliptic  solver,  and  Eq.  (23)  then  yields  the  electron 
velocity  field.  We  then  utilize  exactly  the  same  multidimensional  FCT  transport  algorithm  used 
in  the  previous  section  for  solid  body  rotation  to  integrate  Eq.  (21)  in  time. 

Our  computational  mesh  consists  of  40  grid  points  in  the  x direction,  160  grid  points  in 
the  y direction,  periodic  boundary  conditions  in  both  directions,  and  Ax  = A y = 0.31  km.  Our 
"cloud"  consists  of  a 1-D  gaussian: 

Ne(x,y)  = N0  (1  + lOe'^  ■V°>J/64) 

where  N0  is  the  ambient  background  electron  density  and  y0  is  the  spatial  center  of  the  gaussian 
distribution.  Superimposed  upon  this  distribution  is  a random  x-dependent  perturbation  with  a 
maximum  amplitude  of  3 percent. 

Figures  16-20  show  isodensity  contours  of  Ne/N0  for  the  above  configuration  at  various 
times  in  the  integration.  It  is  seen  that,  as  expected,  the  back  of  the  cloud  (the  upper  half  in 
the  plots)  is  unstable,  growing  linearly  in  the  very  early  stages  of  development.  Non-linear 
effects  soon  enter  the  physics,  however,  as  each  striation  successively  bifurcates,  producing 
smaller  and  smaller  scale  structures,  in  agreement  with  the  results  of  the  ionospheric  barium 
cloud  releases  which  we  are  attempting  to  model.  Two  points  which  bear  on  the  numerics 
should  be  noted:  1)  the  intense  gradients  dictated  by  the  physics  are  nor  diffused  away,  nor  do 
there  appear  in  the  problem  any  of  the  "ripples"  associated  with  numerical  dispersion  which 
normally  appear  when  steep  gradients  try  to  form;  2)  precisely  because  we  did  not  have  to 
resort  to  time-splitting,  none  of  the  usual  time  splitting  phenomena,  such  as  temporal  density 
oscillations  ans  spurious  density  values,  are  evident. 

CONCLUSIONS 

We  have  shown  that  the  algorithm  presented  in  Eq.  (6')  through  (14'),  (17  ) and  (18  ) 
does,  in  fact,  represent  a workable  multidimensional  flux  limiter.  In  addition,  due  to  the  flexi- 
bility in  determining  overshoot  and  undershoot  criteria  inherent  in  the  method,  the  algorithm 
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produces  results  which  are  consistently  equai  or  superior  to  those  produced  using  a time-spli; 


version  of  the  original  flux  limiter  (5),  at  least  for  the  admittedly  limited  class  of  problems 
presented  here. 

For  multidimensional  problems  where  time  splitting  is  unacceptable,  or  for  problems 
where  the  "clipping"  phenomenon  associated  with  the  original  flux  limiter  (5)  is  a serious  prob- 
lem, the  new  algorithm  presented  here  represents  the  only  way  that  FCT  may  be  implemented. 
For  these  problems  the  choice  is  clear,  for  there  is  only  one  option.  Yet  even  in  situations 
where  the  constraints  mentioned  above  do  not  apply,  benefit  may  be  gained  by  implementing 
the  new  algorithm  rather  than  time  splitting  Eq.  (5).  We  do  not  yet  have  enough  experience  to 
give  any  guidelines,  and  can  only  ask  the  prospective  user  to  try  the  method. 

Certainly  the  possiblities  for  modifying  the  basic  scheme  are  endless.  One  could,  for 
instance,  limit  the  antidiffusive  fluxes  only  with  respect  to  maxima,  or  to  minima;  or  he  could 
limit  the  fluxes  sequentially  for  maxima  and  minima,  rather  than  limiting  maxima  and  minima 
simultaneously  in  the  manner  presented  here  (this  last  procedure  will  introduce  an  asymmetry 
between  the  treatment  of  maxima  and  minima  which  may  or  may  not  be  desirable).  Even 
within  time-split  codes  there  are  possibilities.  One  could  time  split  the  one  dimensional  form  of 
the  new  algorithm  rather  than  time  splitting  Eq.  (5);  or  fully  multidimensional  flux  limiting 
could  be  performed  at  the  end  of  each  sweep  of  a time-split  scheme. 

On  NRL’s  Texas  Instruments  ASC  computer,  the  calculations  presented  in  section  VII 
required  93  seconds  and  125  seconds  of  CPU  time  for  the  time-split  and  fully  multidimensional 
cases  respectively,  a cost  penalty  of  slightly  more  than  30%  for  the  multidimensional  limiter. 
Of  course  this  extra  cost  is  highly  problem  dependent.  For  instance  the  striations  code 
described  in  section  VIII  spends  80%  of  its  time  solving  Eq.  (22),  making  the  net  cost  penalty 
of  fully  multidimensional  flux  limiting  only  a few  percent. 


» 


16 


ACKNOWLEDGEMENTS 


1 would  like  to  lhank  Drs.  J.  Boris,  D.  Book  and  J.  Gardner  for  first  introducing  me  to 

FCT,  and  for  innumerable  discussions  over  the  past  few  years  on  the  subject  of  transport  algo- 
rithms. 1 would  also  like  to  especially  thank  Dr.  Book  for  his  indefatigable  encouragement. 

This  work  was  supported  by  the  Defense  Nuclear  Agency. 

REFERENCES 

1.  J.P  Boris  and  D.L.  Book,  Flux-corrected  transport.  I.  SHASTA,  a fluid  transport  algo- 
rithm that  works,  J.  Comp.  Phys.  11  (1973),  38. 

2.  D.L.  Book,  J.P.  Boris,  and  K.  Hain,  Flux-corrected  transport.  II.  Generalizations  of  the 
method,  J.  Comp.  Phys.  18  (1975),  248. 

3 J.P.  Boris  and  D.L.  Book,  Flux-corrected  transport.  III.  Minimal-Error  FCT  algorithms,  J 
Comp.  Phys.  20  (1976),  397. 

4.  D.  Gottlieb,  Strang-type  difference  schemes  for  multidimensional  problems,  SIAM  J. 
Num.  Anal.  9 (1972),  650. 

5.  C.K.  Forester,  Higher  order  monotonic  convective  difference  schemes,  J.  Comp.  Phys.  23 
(1977),  1. 

6.  H.-O.  Kreiss  and  J.  Oliger,  Comparison  of  accurate  methods  for  the  integration  of  hyper- 
bolic equations,  Tellus  24  (1972),  199. 

7.  A.J.  Scannapieco,  S.L.  Ossakow,  S.R.  Goldman,  and  J.M.  Pierre,  Plasma  cloud  late  time 
striation  spectra,  J.  Geophys.  Res.  81(1976),  6037. 


17 


Sixth  order  : 


^ I + ( 1 / 2 ) /A  (/,  + !+/,)  - TT(/,+J  + /,-|) 


+ ^</,+J  + /,-2> 


Eighth  order  : 

'■■*»'» -So 

+ sl/»'+w-s(/-+w 

The  above  fourth  and  eighth  order  forms  are  used  as  the  high  order  (luxes  in  the  main 
body  of  this  paper. 

The  low  order  flux  of  the  leapfrog-trapezoidal  FCT  schemes  is  simply  donor  cell  plus 
a zeroth  order  diffusive  flux  with  coefficient  — . The  donor  cell  algorithm  requires  that 

O 


/ = vw,  where  v is  a convective  velocity.  Specifically, 


, + <]/2)  ~ v / + < 1/2)  wf+l\/2)  g (x/  + i -O  w°)  \i 


v,+d/2)  = y(v,  + v,+1) 

DC  _ 1 < if  v/+(l/2)  > 0 


' + ,,/2>  “ |<+,  if  Vl  + (1/2)  < 0 

0 _ J W"~'  f°r  leaPfr°8  sleP 

W'  |vv,"  for  trapezoidal  step 


A detailed  description  and  analysis  of  these  and  other  high  order  FCT  algorithms  will 
be  discussed  in  a forthcoming  report  by  the  present  author. 
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Fig.  1.  The  eight  possible  configurations  of  the  transported  and 
diffused  solution  w'd  in  the  neighborhood  of  a positive  (rightward- 
directed)  antidiffusive  flux  A,+( U2).  Note  that  configurations  1 through 
4 differ  from  configurations  5 through  8 only  in  the  sign  of  the  quantity 


i-1  i i+1 

x 


Fig.  2.  Perspective  view  of  a two  dimensional  profile  of  the  tran- 
sported and  diffused  solution  w"\  showing  the  four  possible 
anitdiffusive  fluxes  affecting  the  grid  point  (/,  j),  the  directions  of 
which  are  indicated  by  arrows.  Note  that  all  of  the  fluxes  are  paral- 
lel to  the  local  gradient  in  w'd  (as  "antidiffusive"  fluxes  might  be 
expected  to  be),  and  that  w'rf,  is  not  an  extremum.  This  situation  is 
impossible  in  one  dimension,  and  it  is  precisely  this  impossibility 
which  allows  fluxes  to  be  limited  without  regard  to  neighboring 
fluxes  (see  text).  In  two  or  more  spatial  dimensions  a flux-limiting 
formula  must  take  into  account  effects  due  io  multiple  fluxes  acting 
in  concert. 
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Fig.  3.  One  dimensional  profile  of  the  transported  and  diffused 
profile  w'd , showing  the  two  antidiffusive  fluxes  /4,+(l/2)  and 
/4,_( i/j)  whose  collective  effect  must  be  taken  into  account  with 
respect  to  overshoots  and  undershoots  in  the  final  value  of 


Fig.  4.  Comparison  of  old  and  new  flux  limiters  on  narrow  gaussian 
profile  in  passive  convection  for  the  trivial  case  of  a vanishing  velo- 
city field.  The  transport  algorithm  is  LPE  SHASTA.  Note  the 
"clipping"  phenomenon  associated  with  the  old  limiter. 
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Fig.  5.  Same  comparison  as  in  Fig.  4 except  that  the 
velocity  fluid  is  now  finite.  The  profile  has  been  con- 
vected  through  60  grid  points.  Note  the  reduced  clip- 
ping with  the  new  flux  limiter. 


i 

Fig.  6.  Same  comparison  as  in  Fig.  5,  but  with  a more  accurate  tran- 
sport algorithm  (2-8  leapfrog-trapezoidal).  Again  note  the  reduced  clip- 
ping with  the  new  flux  limiter. 
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Fie.  7.  Time  sequence  of  profiles  produced  by  a "perfect  convection 
scheme  acting  on  the  variable  p with  e = 0.2.  The  actual  analytic 
p-ofile  is  shown  as  a solid  line,  and  the  grid  point  values  are  shown  as 
dots.  Note  that  at  time  /„  + 4 At  a grid  point  value  at  (/  + U has 
been  generated  which  is  higher  than  any  grid  point  value  at  the  previ- 
ous time  step.  This  is  the  reason  that  even  the  new  flux  limiter, 
using  Eq.  (17)  and  (18)  for  wmax  and  wmm,  must  still  "clip". 


wtd 

i 


\ 


^ (IF  INTERSECTION 
LIES  BETWEEN  Xj  AND 


/ *i+l) 

\ 


/ 


_i 

i-1 


i + 1 


i+2 


Fig.  8.  A possible  scheme  for  extracting  information  about  extrema 
which  exist  between  grid  points  at  a given  point  in  time.  An  extremum 
is  assumed  to  exist  between  grid  points  i and  i + 1 if  the  intersection  of 
the  right  and  left  sided  extrapolations  of  w"  has  an  x coordinate 
between  x,  and  x,+l.  The  w coordinate  of  the  intersection  is  then  used 
in  the  computation  of  wmax  and  wmm  (see  text). 
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Fig.  9.  Same  as  Fig.  6,  except  that  Eq.  (19)  and  (20),  which  utilize  the 
wpeak  computation  illustrated  in  Fig.  8.,  are  used  to  compute  wmax  and 
wm'n  in  the  new  flux  limiter.  Values  for  the  old  flux  limiter,  since  they 
are  identical  to  those  shown  in  Fig.  6,  are  not  shown.  Note  that  the 
clipping  has  been  virtually  eliminated. 


wtd  PROFILE: 


x 


Fig.  10.  Two  dimensional  profile  of  the  transported  and 
diffused  values  w"\  showing  the  four  antidiffusive 

fluxes  -4, + (1/2).  />  -(1/2).  /.  <4 1,  / + (l/2)>  ancl  -4/.  /— (1/2) 

whose  collective  effect  must  be  taken  into  account  with 
respect  to  overshoots  and  undershoots  in  the  final  value 
of  w"*1-  A perspective  view  of  a similar  profile  is 
shown  in  Fig.  2. 
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-100  CELLS- 


Fig.  11.  Schematic  representation  of  two  dimensional  solid 
body  rotation  problem.  Initially  w inside  the  cut-out 
cylinder  is  3.0,  while  outside  w = 1.0.  The  rotational  speed 
is  such  that  one  full  revolution  is  effected  in  628  cycles. 
The  width  of  the  gap  separating  the  two  halves  of  the 
cylinder,  as  well  as  the  maximum  extent  of  the  "bridge" 
connecting  the  two  halves,  is  5 cells. 


Fig.  12.  Perspective  view  of  initial  conditions  for 
the  two  dimensional  solid  body  rotation  problem. 
Note  that  only  a 50  x 50  portion  of  the  mesh  cen- 
tered on  the  cylinder  is  displayed.  Grid  points 
inside  the  cylinder  have  tv  , = 3.0.  All  others 
have  »f  = 1.0. 
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Fig.  13.  Comparison  of  perspective  views  Fig.  14.  Same  as  Fig.  14,  but  after  628 

of  the  tv  profile  after  157  iterations  (1/4  iterations  (one  full  revolution).  Again  note 

revolution)  with  both  the  old  and  new  flux  decreased  diffusion  with  new  flux  limiter, 

limiters.  The  perspective  view  has  been 
rotated  with  the  cylinder,  so  that  direct 
comparison  with  Fig.  12  can  be  made. 

Again  we  plot  only  the  50  x 50  grid  cen- 
tered on  the  analytic  center  of  the  cylinder. 

Features  to  compare  are  the  filling-in  of  the 
gap,  erosion  of  the  "bridge",  and  the  rela- 
tive sharpness  of  the  profiles  defining  the 
front  surface  of  the  cylinder. 
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PLASMA  DENSITY 
ENHANCEMENT 


b)  t>0 

Fig.  15.  Schematic  representation  of  the  development 
of  a plasma  cloud  (plasma  density  increasing  toward  the 
center)  in  crossed  electric  and  magnetic  fields.  Super- 
imposed on  the  bulk  E„  x B„  motion  is  a steepening  of 
the  rearward  side  of  the  cloud.  This  same  side  is  physi- 
cally unstable  to  small  perturbations. 
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Fig.  16.  Isodensity  contours  of 
plasma  density  at  t = 0 sec.  The  ini- 
tial destribution  for  NJN„  is  a gaus- 
sian  in  y,  centered  at  y — 12.1  km, 
plus  a small  random  perturbation  in  x. 
Contours  are  drawn  for  N1JN„  = 1.5, 
3.5,  5.5,  7.5  and  9.5.  The  area 
between  every  other  contour  line  is 
cross-hatched.  Only  120  of  the  160 
cells  actually  used  in  the  y direction 
are  displayed.  Boundary  conditions 
are  periodic  in  both  directions.  In  our 
plot  B„  is  toward  the  reader,  and  E„  is 
directed  toward  the  right,  and  we 
have  placed  outselves  in  a frame 
moving  with  the  (c/|fl„|2)  E„x  B„ 
velocity.  The  upper  portion  of  the 
gaussian  is  physically  unstable  to  per- 
turbations, while  the  lower  half  is 
(linearly)  stable. 
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T = 3 0 4 SEC 


Fig.  19.  Same  as  Fig.  16,  but  for  t = 
304  sec.  Development  is  fully  non- 
linear, as  the  intense  gradients  and 
associated  high  Fourier  wave  numbers 
become  apparent. 


T = 4 0 7 SEC 


0.0  2.7 

X (KM) 


Fig.  20.  Same  as  Fig.  16,  but  for  t = 
407  sec.  Several  plasma  bifurcations 
are  apparent,  in  agreement  with  the 
experimental  results  from  ionospheric 
barium  cloud  releases,  and  we  have 
maximum  to  minimum  density  varia- 
tions resolved  over  only  2 cells. 
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NASA 

600  Independence  Ave. , S.  W. 
Washington,  D . C.  203^6 
ATTN:  K.  Dub  in’ 


Aerodyne  Research,  Inc. 


Tech/Ops  Building 
20  South  Avenue 


Bui  ling ton 

, MA  01303 

ATTN: 

M.  Camac 

ATTN: 

F.  Bien 

Acroapace 

Corporation 

P.  0.  Box 

92957 

Los  Ange.es,  CA  90009 

ATTN: 

T.  M.  Salmi 

ATTN: 

S . P . Bower 

ATTN: 

V.  Josephson 

ATTN: 

SMFA  for  FWW 

ATTN: 

R.  Grove 

ATTN: 

R.  D.  Rawcliffe 

ATTN: 

T.  Taylor 

ATTN: 

Harris  Mayer 

ATTN: 

D.  C.  Cartwright 

Analytical  Systems  Corporation 
25  Ray  Avenue 
Burlington,  MA  01803 

ATTN:  Radio  Sciences 

Avco-Everett  Research  Laboratory,  Inc. 
2385  Revere  Beach  Parkway 
Everett,  MA  021^9 

ATTN:  Richard  M.  Patrick 

Boeing  Company,  The 
P.  0.  Box  3707 
Seattle,  WA  9812U 

ATTN:  D.  Murray 
ATTN:  Glen  Keister 

Brown  Engineering  Company,  Inc. 
Cummings  Research  Park 
Huntsville,  AL  35807 

ATTN:  David  Lambert  MS  18 

California  at  San  Diego,  Univ.  of 
Building  500  Mather  Campus 
3172  Miramar  Road 
La  Jolla,  CA  92037 

ATTN:  Henry  G.  Booker 
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Calspan 
P.  0.  Box  235 
Buffalo,  N.  Y.  14221 

ATTN:  Romeo  A.  Deliberis 

Computer  Sciences  Corporation 

P.  0.  Box  530 

6565  Arlington  Blvd. 

Falls  Church,  VA  22046 
ATTN:  H.  Blank 
ATTN:  Barbara  F.  Adams 

Comsat  Laboratories 
P.  0.  Box  115 
Clarksburg,  Md.  20734 
ATTN:  R.  R.  Taur 

Cornell  University 

Department  of  Electrical  Engineering 
Ithaca,  N.  Y.  14850 

ATTN:  D.  T.  Farley,  Jr. 

ESL,  Inc. 

495  Java  Drive 
Sunnyvale,  CA  93102 

ATTN:  J.  Roberts 
ATTN:  V.  L.  Mower 
ATTN:  James  Marshall 
ATTN:  R.  K.  Stevens 

General  Electric  Company 
Tempo-Center  for  Advanced  Studies 
8l6  State  Street 
Santa  Barbara,  CA  93102 
ATTN:  Don  Chandler 
ATTN:  DASIAC 
ATTN : Tim  Stephens 


General  Electric  Company 
P.  0.  Box  1122 
Syracuse,  N.  Y.  13201 

ATTN:  F.  A.  Reibert 

General  Research  Corporation 
P.  0.  Box  3587 
Santa  Barbara,  CA  93l°5 
ATTN:  John  Ise,  Jr. 


Geophysical  Institute 
University  of  Alaska 
Fairbanks,  AK  997^1 

ATTN:  Technical  Library 
ATTN:  Neil  Brown 
ATTN:  T.  N.  Davis 

GTE  Sylvania,  Inc . 

I89  B Street 

Needham  Heights,  MA  02194 
ATTN:  Marshall  Cross 

HRB -SINGER,  Inc. 

Science  Park,  Science  Park  Road 

P.  0.  Box  60 

State  College,  PA  16801 

ATTN:  Larry  Feathers 

Honeywell  Incorporated 
Radiation  Center 
2 Forbes  Road 
Lexington,  MA  02175 

ATTN:  W.  Williamson 

Illinois,  University  of 
Department  of  Electrical  Engineering 
Urbana,  IL  6lS0l 

ATTN:  K.  C.  Yeh 

Institute  for  Defense  Analyses 
400  Army-Navy  Drive 
Arlington,  VA  22202 

ATTN:  Ernest  Bauer 
ATTN:  Hans  Wolfhard 
ATTN:  J.  M.  Aein 
ATTN:  Joel  Bengston 

Inti  Tel  & Telegraph  Corporation 
500  Washington  Avenue 
Nutley,  N.  J.  07H0 

ATTN:  Technical  Library 

ITT  Electro-Physics  Laboratories,  Inc. 
91i0  Old  Annapolis  Road 
Columbus,  Md.  21043 

ATTN:  John  M.  Kelso 
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Johns  Hopkins  University- 
Applied  Physics  Laboratory 
6621  Georgia  Avenue 
Silver  Spring,  MD  20910 

ATTN : Document  Librarian 

Lockheed  Missiles  & Space  Co.,  Inc. 

P.  0.  Box  504 
Sunnyvale,  CA  9^38 
ATTN:  Dept.  60-12 

Lockheed  Missiles  and  Space  Company 
3251  Hanover  Street 
Palo  Alto,  CA  94304 

ATTN:  Billy  M.  McCormac,  Dept  52-14 
ATTN:  Martin  Walt,  Dept  52-10 
ATTN:  Richard  G.  Johnson,  Dept  52-12 
ATTN:  JOHN  CLADIS 


MIT  Lincoln  Laboratory 
P.  0.  Box  73 
Lexington,  MA  02173 

ATTN:  Mr.  Walden,  X113 

ATTN:  D.  Clark 

ATTN:  James  H.  Pannell,  L-246 

ATTN:  Lib  A-082  for  David  M.  Towle 

Martin  Marietta  Corporation 
Denver  Distribution 
P.  0.  Box  179 
Denver,  CO  80201 

ATTN:  Special  Projects  Program  248 

Maxwell  Laboratories,  Inc. 

9244  Balboa  Avenue 
San  Diego,  CA  92123 

ATTN:  A.  J.  Shannon 
ATTN:  V.  Fargo 
ATTN:  A.  N.  Rostocker 

McDonnell  Douglas  Corporation 
5301  Bolsa  Avenue 
Huntington  Beach,  CA  92c57 
ATTN:  J.  Moule 
ATTN:  N.  Harris 
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Mission  Research  Corporation 
735  State  Street 
Santa  Barbara,  CA  93101 
ATTN:  R.'  Hendrick 
ATTN:  Conrad  L.  Longmire 
ATTN:  Ralph  Kilb 
ATTN:  R.  E.  Rosenthal 
ATTN:  R.  Bogusch 
ATTN : David  Sowle 
ATTN:  M.  Scheibe 
ATTN:  P.  Fischer 
Mitre  Corporation,  The 
Route  62  and  Middlesex  Turnpike 
P.  0.  Box  208 
Bedford,  MA  01730 

ATTN:  Chief  Scientist  V/.  Sen 
ATTN:  S.  A.  Morin  M/S 
ATTN:  C.  Ilirding 


North  Carolina  State  Univ  At  Raleigh 
Raleigh,  N.  C.  27507 

ATTN:  SEC  Officer  for  Walter  A.  Flood 

Pacific-Sierra  Research  Corp. 

1456  Cloverfield  Blvd. 

Santa  Monica,  CA  90404 

ATTN:  E.  C.  Field,  Jr. 

Philco-Ford  Corporation 
Western  Development  Laboratories  Div 
3939  Fabian  Way 
Palo  Alto,  CA  94303 

ATTN:  J.  T.  Mattingley  MS  X22 

Fnotometrics,  Inc. 

442  Marrett  Road 
Lexington,  MA  02173 

ATTN:  Irving  J.  Kofsky 

Mitre  Corporation,  The 
Westgate  Research  Park 
1820  Dolley  Madison  Blvd. 

McLean,  VA  22101 

ATTN : Allen  Schneider 
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Physical  Dynamics,  Inc. 

P.  0.  Box  1069 
Berkeley,  CA  9**701 

ATTN:  Joseph  B.  Workman 

Physical  Sciences,  Inc. 

607  North  Avenue,  Door  18 
Wakefield,  MA  01880 
ATTN:  Kurt  Wray 

R & D Associates 

P.  0.  Box  3580 

Santa  Monica,  CA  ^ObO^ 

ATTN:  Robert  E.  Lelevier 
ATTN:  Forest  Gilmore 
ATTN:  Richard  Latter 
ATTN.  William  B.  Wright,  Jr. 

R & D Associates 
1815  N.  Ft.  Myer  Drive 
11th  Floor 
Arlington,  VA  2209 

ATTN:  Herbert  J.  Mitchell 

Rand  Corporation,  The 
1700  Main  Street 
Santa  Monica,  CA  9 0*406 

ATTN:  Cullen  Crain 


Science  Applications,  Inc. 

P.  0.  Box  2351 
La  Jolla,  CA  92038 

ATTN:  Daniel  A.  Hamlin 
ATTN:  D.  Sachs 
ATTN:  E.  A.  Straker 

Space  Data  Corporation 
1331  South  26th  Street 
Phoenix,  AZ  8503^ 

ATTN:  Edward  F.  Allen 
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Stanford  Research  Institute 
333  Ravenswood  Avenue 
Menlo  Park,  CA  9^025 
ATTN:  M.  Earon 
ATTN:  L.  L.  Cobb 
ATTN:  Walter  G.  Chestnut 
ATTN:  David  A.  Johnson 
ATTN:  Charles  L.  Rino 
ATTN:  E.  J.  Fremouw 
ATTN:  Ray  L.  Leadabrand 
ATTN:  Donald  Neilson 

Stanford  Research  Institute 
306  Wynn  Drive,  N.  W. 

Huntsville,  AL  35805 

ATTN:  Dale  H.  Davis 

Technology  International  Corporation 
75  Wiggins  Avenue 
Bedford,  MA  01730 

ATTN:  W.  P.  Boquist 

TRW  Systems  Group 
One  Space  Park 
Redondo  Beach,  CA  90278 
ATTN:  P.  H.  Katsos 
ATTN:  J.  W.  Lowery 

Utah  State  University 
Logan,  UT  8U321 

ATTN:  C.  Wyatt 
ATTN:  D.  Burt 
ATTN:  Kay  Baker 
ATTN:  Doran  Baker 

Visidyne,  Inc. 

19  Third  Avenue 

North  West  Industrial  Park 

Burlington,  MA  01803 

ATTN:  William  Reidy 
ATTN:  Oscar  Manley 
ATTN:  J.  W.  Carpenter 
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